function output = rhythm_gpu(y)
    if nargin < 1, y = wavread('technobeat.wav'); end
       
    bandlimits=[0 200 400 800 1600 3200];
    maxfreq=44100;
    y = gsingle(y);
    
    yFFTed = fft(y);
    
    n = length(yFFTed);
    nbands = length(bandlimits);
    
    for i = 1:nbands-1
        bl(i) = floor(bandlimits(i)/maxfreq*n/2)+1;
        br(i) = floor(bandlimits(i+1)/maxfreq*n/2);
    end
    
    bl(nbands) = floor(bandlimits(nbands)/maxfreq*n/2)+1;
    br(nbands) = floor(n/2);
    
    output = gzeros(n,nbands);
    
    
    
    for i = 1:nbands
        output(bl(i):br(i),i) = yFFTed(bl(i):br(i));
        output(n+1-br(i):n+1-bl(i),i) = yFFTed(n+1-br(i):n+1-bl(i));
    end
    
    output(1,1)=0;
    
    
    
    winlength = .4;
    
    hannlen = winlength*2*maxfreq;
    
    hann = [zeros(n,1)];
    
    
    for a = 1:hannlen
        hann(a) = (cos(a*pi/hannlen/2)).^2;
    end
    
    
    wave = gones(size(output));
    gfor k = 1:nbands
    wave(:,k) = real(ifft(output(:,k)));
    gend
    freq = gones(size(wave));
    
    
    wave = single(wave);
    for i = 1:nbands
        for j = 1:n
            
            if wave(j,i) < 0
                wave(j,i) = -wave(j,i);
            end
        end
        tempwave = gsingle(wave(:,i));
        freq(:,i) = fft(tempwave);
        
    end
    
    
    hann = gsingle(hann);
    freq = gsingle(freq);
    filtered = gzeros(size(output));
    output = gzeros(size(output));
    gfor m = 1:nbands
        filtered(:,m) = freq(:,m).*fft(hann);
        output(:,m) = real(ifft(filtered(:,m)));
    gend
